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A  FAST  PARABOLIC  MODULE  FOR  THE  SOLUTION  OF  MHD  CHANNEL  FLOW  EQUATIONS 

BETWEEN  ELECTRODE  WALLS 


by  J.P.F.  Lindhout 
H.  Snel 
W.F.H.  Merck 


INTRODUCTION 

The  present  study  is  part  of  the  NLR  research  program  aimed  at  generating 
numerical  solutions  of  Magnetic  Hydrodynamic  (MHD)  duct  flow.  This  research 
program  is  incorporated  in  the  investigations  into  plasma-flow  phenomena, 
which  are  performed  in  the  Group  Direct  Energy  Conversion  of  the  EUT. 

MHD  energy  conversion  is  a  process  in  which  the  motion  and  the  thermal 
energy  of  an  electrically  conducting  fluid  through  an  applied  magnetic 
field  perpendicular  to  the  motion  is  converted  into  electrical  energy. 

This  process  yields  higher  conversion  efficiences  than  conventional 
processes.  A  schematic  view  of  an  MHD  generator  channel  is  given  in  fig.  1. 
The  channel  is  horizontally  bounded  by  insulator  walls  and  in  the  vertical 
direction  by  electrode  walls  in  which  segmented  electrodes  are  placed. 

The  magnetic  field  B  is  applied  perpendicular  to  the  insulator  walls.  The 
induced  electrical  field  E  equals  the  vector  product  uxB.  It  produces  a 
current  in  the  conducting  medium  which  flows  through  the  electrodes. 

The  complete  system  of  equations  describing  the  phenomena  can  be  divided 
into  three  groups: 

i)  the  gasdynamical  -  or  global  equations 

ii)  the  electron— gas  equations 

iii)  the  electromagnetic  field  equations. 

The  first  two  groups  of  equations  are  of  parabolic  type.  Parabolic  equa¬ 
tions  in  fluid  dynamics  frequently  characterize  situations  in  which  the 
dependent  variables  have  large  gradients  normal  to  the  main  flow  direction. 
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For  this  type  of  equations,  two  fourth-order  accurate  solution  methods 
are  presented.  Both  methods  are  based  on  cubic-spline  representation  of 
the  dependent  variables  normal  to  the  marching  direction.  A  stability 
analysis  will  be  given  and,  to  show  the  accuracy,  both  methods  are  applied 
to  a  simple  flow  problem,  of  which  the  exact  solution  is  known. 

One  of  the  fourth-order-spline  methods  is  used  to  solve  the  electron-gas 
equations  which  have  substantially  more  complex  boundary  conditions  than 
the  global -gas  equations.  Spline  methods  are  preferable  in  situations  like 
this  where  derivatives  enter  the  boundary  conditions. 

A  second-order  difference  scheme  was  already  applied  by  the  present  authors 
to  solve  the  global-gas  equations  between  insulator  walls  (ref.  l). 

2.  CUBIC-SPLINE  SOLUTION  OF  PARABOLIC  DIFFERENTIAL  EQUATIONS 

Rubin  and  Graves  (ref.  2)  introduced  the  application  of  collocational  cubic 
splines  to  the  solution  of  parabolic  partial  differential  equations.  For 
two-dimensional  problems,  conventional  finite  differencing  took  place  in  * 
the  marching  direction.  The  resulting  two— point-boundary-value  problem,  nor¬ 
mal  to  the  marching  direction,  was  solved  by  cubic  splines.  The  spline  is 
required  to  satisfy  the  second-order  differential  equation  at  each  grid 
point.  Applying  the  spline  continuity  relations  one  can  obtain  a  block  tn- 
diagonal  system  of  equations  for  the  function  values,  say  u,  its  first 
(m)  and  second  derivatives  (M)  in  the  normal  direction.  This  method  is  of 
second-order  accuracy  in  u  and  m,  because  errors  in  M  itself  are,  according 
to  the  cubic— spline  approximation,  proportional  to  the  square  of  the  inter¬ 
val  width. 

Daniel  and  Swartz  (ref.  3)  suggested  a  method  of  obtaining  fourth-order 
accuracy  for  two-point-boundary-value  problems,  basically  by  eliminating 
the  second-order  truncation  error  in  M,  which  is  h“(cT*u/3y  )/l2+0(h  ).  For 
the  inner  points,  the  leading  term  of  the  truncation  error  was  approximated 
by  means  of  second-order  central  differences  of  M.  At  the  boundary  points 
non  central  differences  had  to  be  used.  The  collocation  equations  contained 
information  of  other  points  and  the  resulting  equations  could  only  be  re¬ 
duced  to  a  five— diagonal  system. 

Rubin  and  Khosla  (ref.  4)  explored  this  method  extensively  for  parabolic 
equations. 
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To  demonstrate  the  details  of  both  proposed  methods,  the  boundary  layer 
equations  of  an  incompressible  fluid  near  a  stagnation  point  are  solved. 

The  governing  flow  equations  can  be  written  as  follows: 

2  2 

conservation  of  momentum:  u(Du/Dx)+v(<)u/£)y)=ue(due/dx)+V(3  u/3y  ),  (l) 

conservation  of  mass  :  (<)u/c)x)  +  (5v/0y)  =0,  (2) 

in  which  x  and  y  form  a  Cartesian  coordinate  system,  u  and  v  are  velocities 
in  x-  and  y-direction,^  is  the  known  viscosity,  ug  is  the  given  velocity  at 
the  edge  of  the  boundary  layer.  When  equation  (l)  is  used,  equation  (2)  can 
be  replaced  by: 

?)  (v/u)/fcy=  -  [ue(due/dx)+^(02u/()y2)]/u2,  (3) 

in  which  only  y— derivatives  of  the  unknowns  occur.  The  boundary  conditions 
are : 


y  =  0  :  u  =  0,  v  =  0  (4a) 

y  -co  :  u  =  ug  (4b) 

When  stagnation  conditions  are  assumed,  i.e.  (due/dr)  is  constant,  a  similar¬ 
ity  solution  of  (l)  and  (2)  exists,  which  permits  the  determination  of  the 
accuracy  of  the  proposed  methods  (ref.  5)  •  A  grid  is  assumed  at  which  a 
grid  function  u^  exists:  u?»u(x=i.k,  y=(j-l).h).  Equation  (l)  is  centered  at 
the  grid  point  (i+1,  j).  The  derivative  in  the  marching  direction  is  approx¬ 
imated  by  the  second— order  difference  formula: 

•0u/j>x->,(3u^+1  -  4u^  +  u?1_1)/(2k),  ^ 

which  is  substituted  in  (l).  The  non-linear  terms  at  the  level  i+1,  in  the 
resulting  equation  are  linearised:  (u^+1)2  by  the  Newton-Raphson  method  and 
for  in  v^+^  (Ou/Oy)  the  value  of  the  last  iteration  is  taken. 

Equation  (l)  can  be  written  with  the  help  of  equation  (5)  and  the  linearizing 
techniques,  as  an  ordinary  differential  equation: 

a^  h2(02u/®y2)^+1  +  b  Jitou/cy)^.  +  c..  u^  =  djf  (6) 

(i-1,2,...,  j=1 ( 1 )N ,  N  odd) 

in  which  a,  b,  c,  d  are  functions  of  h,  U(>  (du^/dx)  and  the  last  iterate  of 
the  unknowns  u  and  v. 
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At  each  collocation  point  the  following  approximations  are  made: 

(3u/ay)^+1~m^+1  ,  (g2u/ay2)^+1^M^+1,  (7a, b) 

where  m  and  M  are  the  first  and  second  derivative,  resp.  of  the  cubic  spline. 
According  to  Ahlberg,  Nilson  and  Walsh  (ref.  6),  the  first  approximation  is 
fourth-order  accurate  and  the  last  one  second-order  accurate. 

From  the  equidistant  cubic-spline  approximation  of  u,  it  is  possible  to  de¬ 
duce  the  following  four  indepent  equations  for  u,  m  and  M  at  three  adjacent 


grid  points  (ref.  6): 

M'5-1  =  6  (u^-u^-1  )/h2-2(m^+2m^-1  )/h,  (8a) 

=  3(u^+1-2u^+u^  ^  )/h_-(m''+'''-m^  (8b) 

=  -6(u^+1-u^)/h2+2(2m^+1-Hn^)/h  (8c) 

uj+1-uj  1  =  h(m'’+d+4m^+m^  ^)/3»  (8d) 


for  convenience  the  subscript  i+1  is  omitted. 


These  equations  are  employed  to  approximate  m  and  M.  To  demonstrate  this 
we  write  equation  (6)  for  j-1 ,  j  and  j+1  together  with  equation  (8)  as  a 
vector  equation: 

A/1  +  B  UJ  +  C  UJ+1  =  D  ,  j  =  2  (1)  N-l  (9) 

J  J  J  J 

in  which:  =  (h2M^,  hm^,  u^)T 

Dj  =  (dj-l’  dj’  dj+l’°’  °’  °»  °)T’ 


and 


By  prenmlti plication  of  (9)  by  a  suitable  row  vector  r  =  (r^,..r„)  it  is 
possible  to  obtain  a  scalar  equation  with  only  three  unknowns.  In  order 
to  extract  from  (9)  on  equation  in  which  only  the  function  values  u^  , 
uj  and  u^+1  occur  and  no  derivatives,  r  has  to  satisfy: 

1  O  I  O  1  O 

rA .  =  rA  =  rB  =  rB  .  =  rC  =  rC  =  0  ( 

3  3  3  3  3  3 

in  which  the  superscript  the  involved  column  indicates. 

Bow  r  is  fully  determined  hy  (10),  except  for  a  multiplicative  constant, 
and  after  multiplication  of  (9)  by  the  chosen  r  we  obtain: 


a .  u^  1  +  |3  .  u^  +  y .  u3+1  =  p  .  j  =  2(1  )N 

J  J  J  J 


(Ha) 


Mixed  boundary  conditions  of  the  type: 


ah  M1+bhm,+cu-1=d 
o  1  o  1  o  1  o 


can  be  handled  by  the  same  technique.  System  (9)  is  extended  with  the  boun¬ 
dary  condition  by  extending  the  matrices  A,  B,  C  and  vector  D  by  appropriate 
elements.  By  premultiplication  of  a  new  row  vector  r,  now  with  eight  elements, 


one  can  obtain: 


Plu  +  Y1U  =  P1 


r  has  to  satisfy  (10)  and  moreover: 


rC  .  =  0. 
3 


(lib) 


Mixed  boundary  conditions  at  the  upper  boundary  can  be  treated  equivalently 


and  lead  to: 


*  v" 


(11c) 


The  tri -diagonal  system  (ll )  can  be  solved  with  standard  algorithm’s  for  u^. 
The  m^  can  be  solved  from  equation  (6)  and  (8c).  An  initial  condition  for 
m  at  the  lower  boundary  is  obtained  in  the  same  fashion  as  (lib).  M  in  its 
turn  follows  from  the  original  differental  equation.  All  quantities  are 
derived  second  order  accurate  in  h. 


To  obtain  fourth-order  accurate  results  equation  (6)  is  solved  once  again 
with  stepsize  2h  in  the  y-direction,  while  only  the  odd  points  of  the  grid 
are  used. 

The  first  solution  is  indicated  by  *  and  by-v>  the  solution  in  which  only 
the  odd  points  are  used.  The  leading  term  of  the  truncation  error  of  the 
results  of  odd  points  is  four  times  as  large  as  the  results  of  the  first 


solution  for  all  points.  Richardson-extrapolation  consists  in  finding  a 
more  accurate  result  at  the  odd  points  by  elimination  of  the  leading  trun¬ 
cation  term: 

P-(4f*-p/3,  (13) 

in  which  f  stands  for  u^+1 ,  (Ou/0y)^+1  ,  and  (32u/E)y2)^+1  ,  j=l(2)N. 

In  the  further  computations,  it  is  necessary  to  obtain  results  at  the  even 
points  as  accurate  as  those  at  the  odd  points  along  a  normal.  This  is  accom¬ 
plished  in  the  following  way: 

The  function  values  and  the  first  and  second  derivatives  at  the  odd  points 
hold  exactly  enough  information  to  define  a  quintic  hermite  polynomial  over 
the  intervals  between  each  pair  of  successive  odd  points. 

At  the  even  points,  the  quintic  polynomial  is  used  to  interpolate  the  un¬ 
knowns  : 

uj=(uJ-1+u'i+1  )/2+5h[(0u/Sy)  j-1-(t>u/0y)  j+1]/l6+h2{(d2u/c>y2)  j+1  + 

(c)2u/Sy2)J-1}/l6,  j=2(2)N-1,  ( 14) 

(Ou/®y)  ‘U— 15(u'j-1-uj+1  )/l6h-7{(«>u/c>y)  j_1  +  (t>u/Dy)  ^+1J/l6+ 

-h{(S2u/0y2)J_1-O2u/3y2)J+1]/l6,  j=2(2)N-1,  (15) 

(O2u/0y2)  j=-3((t>u/Dy)  J"  -(Ou/0y)  j+1J/4h-((®2u/&y2) j-1+ 

+  (t)2u/Dy2)j+1}/4,  j=2(2)N— 1.  (16) 

The  idea  of  applying  Richardson-extrapolation  to  improve  the  results  of 
numerical  techniques  is  not  new.  Richardson-extrapolation  is  employed  for 
several  years,  e.g. ,  in  the  box  method  of  Keller  (ref.  7),  for  the  compu¬ 
tation  of  2D  and  3D  boundary  layer  flows. 

The  new  feature  here  is  that,  in  the  case  of  the  cubic-spline  solution,  it 
is  permitted  to  apply  Richardson-extrapolation  for  each  marching  step  before 
the  computation  is  proceeded  to  the  next  x-station. 

The  second  method  presented  here  eliminates  the  truncation  term  h  (a  u/oy*+)/l2 

of  equation  (7b).  By  differentiating  the  original  equation  (l)  twice, 

4/4 

we  can  express  0  u/Oy  in  second  and  lower  derivatives,  and  we  can  appro¬ 
ximate  the  second  derivative  in  (l)  at  each  collocation  point  by: 
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tfVcrV'v  M-h2  {(uv/v»+m)  (t>rn/Dx)  +  (viVv)-?)u/c>x)M+u('SM/c>x)}/l  2>>+0(h4 )  ( 17) 

The  x-derivatives  in  (17)  are  approximated  in  the  same  way  as  described  al¬ 
ready  for3u/bx.  The  non-linear  terms  are  linearized  in  the  Newton-Raphson 
fashion,  except  for  the  terms  in  which  v  appears,  which  are  taken  from  the 
most  recent  iteration.  The  resulting  linear  equation  in  u,  m,  and  M  can 
further  be  treated  in  exactly  the  same  way  as  in  the  cubic-spline  case.  It 
is  necessary  to  store  also  m  and  M  at  three  successive  x- stations.  Unlike  the 
method  of  ref.  (3)  the  use  of  unsymmetrical  differences  at  the  wall  is  avoided; 
this  can  be  an  important  advantage  in  the  case  of  steep  gradients.  The  fact 
that  a  tri-diagonal  system  of  equations  results  from  this  method  is  also  a 
clear  advantage.  However,  the  differentiation  can  be  a  tedious  task.  For  com¬ 
plicated  equations  the  use  of  computerised  formula  manipulation  can  be  of 
great  help. 

At  this  point,  it  is  clear  that  either  of  the  two  methods  have  yielded  values 
for  u  and  its  derivatives.  The  RHS  of  equation  (3)  can  be  computed  and 
equation  (3)  can  be  written  as  a  first-brder  differential  equation  in  the 
unknown  (v/u).  The  solution  procedure  can  be  the  same  as  that  used  for  the 
cubic— spline  solution  of  (6)  if  an  additional  boundary  condition  is  given. 

It  is  consistent  with  the  asymptotic  behaviour  of  (l)  and  (2)  to  apply: 

y  =  (N-I).h  :  ZT  (v/u)/c)y2  =0  (18) 

Because  in  (3)  the  second  derivative  is  missing,  the  solution  of  (3)  is 
fourth— order  accurate. 

For  several  stepsizes  h  in  y-direction,  computations  are  performed  starting 
at  x=0.1  untill  x=1.0,  k=0.01.  The  error  exhibited  at  the  last  x-station  in 
the  friction  coefficient  c^  is  depicted  in  fig..  2  for  the  two  methods  present¬ 
ed.  For  comparison  results  obtained  by  apli cation  of  the  cubic-spline  method 
and  a  simple  difference  scheme  are  included.  From  the  slope  of  the  curves, 
it  can  be  inferred  that  the  cubic-spline  method  and  the  difference  method 
are  second-order  accurate  indeed.  Both  methods  proposed  are  fourth-order 
accurate,  except  for  the  first  one  in  the  case  of  a  very  small  number  of 
netpoints. 

3.  STABILITY  ANALYSIS 


We  examine  the  stability  properties  of  the  presented  fourth-order  methods 
using  the  simple  parabolic  equation: 


u  =<y  u 
x  yy 


(19). 
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3  loes  not  imply  a  severe  restriction  of  the  analysis,  because  the  highest 
Lvatives  are  retained. 

i  is  discretized  according  to  (6),  the  following  equation  at  each  collo- 
Lon  point  is  obtained: 


3u^+r2k *(*  u/fty" ) i+ji=4u^-u^_1  • 


3  equation  can  be  solved  by: 

i)  the  cubic— spline  method,  for  short:  CS— method 

ii)  the  method  in  which  two  cubic-spline  solutions  are  used  to  re¬ 
duce  the  truncation  error  by  Richards on-extrapolation,  followed 
by  quintic  hermite  interpolation  for  the  even  points,  for  short: 
RH-method 

iii)  the  method  in  which  the  truncation  error  is  reduced  by  estimation 
of  the  truncation  error  by  differentiation,  for  short:  D-method. 


13-method :  Application  of  equations  (7b )-(8c) permits  the  elimination  of 
n's  and  M’s,  which  gives  rise  to  the  following  equation  in  u: 
if,,  j- h.,  ,  J  -i-n  - 


4(u4  l4'2+u^+1  )-(u,J~'+4uJ  +u^TI) 

'  1  1  1  '  '  1—1  1-1  1-1 


?he  stability  of  this  scheme  can  be  investigated  by  the  von  Neumann 
lethod.  Because  (21)  is  a  three-level  scheme,  the  analysis  leads  to  a 
;econd— order  amplification  matrix  (ref.  8).  The  Pourier-decomposition  of 

i  can  be  written  as: 

u?  =  £  1  exp  (ijnh)  (22) 

n  which  i  is  the  imaginary  unit  and  n  the  wave  number.  Substitution  of 
22)  into  (21)  gives  a  scalar  equation  in  £l  ^ ,  £l^  and  £l^  ^  which 
an  bo  transformed  into  a  vector  equation  by  the  introduction  of  ^>2: 


§2.  .  =  £l. 

J  1+1  i 


’he  resulting  vector  equation  is  given  by: 

i  i+1  =  A  li  »  . 

4/r  -1/1)  j  £  A  T 

amplification  matrix  A  =  I  I  and  3  =(  J2)  , 

\  1  0  / 


. /,  /.  2v«-/ . 


..  \  //«  Ui  \ 
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The  eigenvalues  of  A  satisfy  the  following  quadratic  equation: 

rX'  —  4  A  +  1  =  0  (25) 

If  3^r<4,  the  roots  of  (25)  are  positive  and  smaller  than  or  equal  to 
one.  If  r^4,  the  roots  are  the  complex  conjugate  of  each  other,  and  it 
follows  at  once  that: 

lx,|2 -U2|2  -  A,  X,  -\A2-7r<74  (*«> 

Consequently  it  can  he  concluded  that  the  CS-method  is  unconditionally 
stable. 

ii)  The  RH-method:  The  results  of  the  CS— method  at  all  points,  with  step- 
size  h,  are  indicated  by  a  and  the  solution  at  odd  points,  stepsize  2h  is 
indicated  by»v  •  The  amplification  matrix  for  the  fourth-order  accurate 
results  at  the  odd  points  is  constructed  according  to  equation  ( 1 3) : 

4(4/r*  -  7?)  -73(4/r*  -  7?) 

1  0 

Its  eigenvalues  satisfy  also  the  quadratic  equation  (25) ;  the  leading 
term  can  then  be  expressed  as: 

(28) 

With  the  same  reasoning  as  above  we  can  proof  that  the  moduli  of  the 
eigenvalues  are  equal  to  or  less  than  one. 


(k/h  )g  (1-cos1 


(2+cos ’ 


odd  1+5(k/h2 )<y(  1-cos  f)/( 2+cos  f  ) 


To  examine  the  stability  of  the  even  points,  a  Fourier  decomposition  of 
m  and  M  of  the  cubic-spline  solution  is  assumed: 


mi  =  71  i  exP(IJ)(')/h» 
Mi  =  ^i  ^P^j/O/h2 


(29) 


From  equations  (8),  the  following  relations  can  be  derived  between 
the  Fourier-coefficients  of  m  and  M  and  those  of  u: 


W  =  3Isin /(2+cos  ^  )  €l, 

!?1  =  -6(l-cos  f)/(2+cosf  )  4l 


(30) 
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It  can  be  shown  that  these  relations  also  hold  for  the  fourth— order  ac¬ 
curate  results  at  the  odd  points.  Substitution  of  the  three  Pourier-de— 
compositions  and  equation  (30)  into  equation  (14)  gives  the  connection 
between  the  eigenvalues  of  the  odd  and  the  even  amplification  matrix: 

l>wj  =  t5-(5-o°s  f  )2  /  8}/(2+cos/')  |Aodd|  (31) 

Hence  the  moduli  of  the  even  eigenvalues  are  equal  to  or  less  than  the 
odd  ones. 

In  fig.  3»  the  eigenvalues  are  depicted  for  the  three  cases  treated  in 
this  chapter.  The  conclusion  can  be  drawn  that  the  RH-method  is  uncon¬ 
ditionally  stable. 

2  2 

iii)  D -met hod:  The  introduction  of  a  better  approximation  for  t>  u/py  with 
respect  to  the  CS— method  does  not  change  the  stability  analysis  because 
only  higher-order  terms  are  added,  which  have  no  influence  upon  the 

stability. 


4.  THE  ELECTRON  GAS  EQUATIONS 

The  RH-method  is  used  to  solve  the  electron  gas  equations  along  the  electrode 
wall  of  an  MHD  channel.  This  exercise  seems  to  be  of  sufficient  difficulty 
to  test  the  method  in  practical  circumstances. 

The  electron  gas  equations  are  derived  from  the  physical  model  of  and  MHD- 
duct,  as  used  for  electric  power  generation,  with  argon  as  a  working  medium 
seeded  with  a  small  fraction  of  cesium  vapour  to  increase  the  electrical 
conductivity .  The  current-carrying  electrons  are  created  by  ionization  of  the  Cs , 
consequently  the  plasma  consists  of  electrons,  CS-ions ,Cs-atoms  and  Ar-atoms.  Be¬ 
cause  there  is  little  interaction  between  the  electrons  and  the  abundant  Ar- 
atoms,  the  electrons  behave  as  a  separate  gas,  with  specific  electron  gas 
temperature  T^  (ref.  9j  10,  11,  12). 

The  interaction  of  the  Ar-Cs  plasma  with  the  applied  magnetic  field  introduces 
a  Lorentz  force  in  the  global  momentum  equation,  and  enthalpy  extraction  and 
ohmic  dissipation  terms  in  the  global  energy  equation.  These  equations  are 
not  dealt  with  in  this  paper;  the  gasdynamic  quantities  p,  u,  v,  T  are  sup¬ 
posed  to  be  known. 


Here  the  solution  of  the  electron  gas  equations  is  emphasised,  because  they 
offer  the  heaviest  problems  in  numerical  computations.  These  problems  are 
generated  by  the  steep  gradients  at  the  wall  and  by  the  exponential  relations 

between  the  coefficients  of  the  differential  equations  and  the  electron 
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Morecver  the  electron  gas  tends  to  be  unstable,  which  is  caused  by  a  positive 
feedback  process  called  ionization  instability.  This  can  be  dealt  with  by  the 
use  of  effective,  time-averaged,  values  for  electrical  conductivity  (3'rf.(  ) 
and  Hall  parameter  (/3eff)  (ref.  11). 

In  our  model,  turbulent  terms  are  omitted  and  the  electrode  wall  is  infinitely 
segmented,  which  results  in  an  axial  current  component  Jx=0.  Non-equilibrium 
ionization  and  ambipolar  diffusion  are  taken  into  account  in  the  electron  gas 
equations,  which  are  merely  stated  here  (ref.  13,  14) J 

continuity 

u(DCe/Ox)  +  v  QC^y)  =  meSe  4  Jp/s)  {  Cjl+T/T) }  (32) 


energy 


/°  cv,e  Ce[^Te/Qx)+v(dTe/dy)J  -  J  - 

Jy/°reff"Hb(Te-T)  +  (®^y)  [  3(k/e)2TeG'eff(^Te/by)+2.5kTe/^/(meS), 
T)  (  Ce(1+Te/T)]  /t>yj  -1.5kTeO/£>y)[y^/(meS)o[ce(l+Te/T)}^yJ  + 


-  n  (Et+c  m  T  ) 
ex  I  v,e  e  e ' 


(33) 


The  symbols  used  in  these  equations  are: 

Cg  =  n^m^^o  ,  electron  gas-concentration 
n^  =  number  density  of  electrons 

=  ne^ns_ne^6,22~n8  exp (-2. 556  e/kTj  + 


'10 


-11^2.58^^  exp(l  .337e/kTg) ,  electric  source  term 


m  =  mass  of  an  electron 
e 

n  =  number  density  of  seed  particles 
s 

n  =  number  density  of  argon 

3. 

u,  =  viscosity 

Ej  =  ionization  energy  of  Cs 

e  =  charge  of  an  electron 

k  =  constant  of  Boltzmann 
2 

J ~  ohmic  dissipation 
J  =  current  density 


A, .(-1-5,6'-  y 

eKc 


otherwise  G”eff  =  =  eB/m^ 
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~  Sr  ^§e^1+Te//T^I  =  ambi-P°lar  ion  flux 

rt  ^r>  ,  f  ^e,i+*e,as  v^efaa|  .  . 

H,  =  3/0 C  k  — 2 - 2 —  +  — 2 —  ,  collision  term 

d  /  e  L  m  rn  J  1 

s  a 

m_  =  mass  of  an  argon  atom 
a 

m„  =  mass  of  a  seed  particle 

O 

S  »  Schmidt  number 

\>  =  V1  +  V>  +  V 

c  e,i  K e  ,as  e,aa 

\>  ■  =  1 • 55inn  T  2  In l  8.76  ^(T  ^/n  )eJ  ,  collision  frequency  electron- 

0  1 1  iuee  (  1  u  6  e) 

Cs-ion 

—20  ~ 

V  n  (8kT  /tt  m  )s,  collision  frequency  electron-Ar-atom 

e  1  3,0,  I  u  a  e  e 

—17  — 

V>  =  i*n  (n  -n  )(8kT  /rr m  )3,  collision  frequency  electron-Cs-atom 
e  1  a  s  1  u  se  s'  e 

The  energy  transport  term  perpendicular  to  the  wall  contains  electron  heat 
conduction,  heat  transport  by  electric  current  and  ambi-polar  diffusion. 


The  boundary  conditions  are  at  the  center  line  of  the  duct : 

Ocjdy  =  0Te/9y  =  0 


(34) 


At  the  wall,  a  two-layer  model  is  used.  It  is  assumed  that  adjacent  to  the 
wall  a  collision— less  electrostatic  sheath  exists  which  is  directly  bounded 
by  the  continuum  plasma  (ref.  13,  1 4 »  15»  16).  This  model  leads  to  the  wall 
boundary  conditions: 

(T+Te)(0Ce/0y)  +  ((c)Te/3y)  -  -^(<>T/Dy)-/£g  (^"P  t]  Cg  =  0  (35) 

and: 


fVJ*  A)I^[P  ^  +  3  5T  ■  0 


for  y= 0,  where  the  sheath  voltage  drop 


kT  r  Jv  J 

Af  - - 2.  in  [  .LJ 

•  e  L  n  e 


Tm 


T  m. ' 
e  l 


k] 


(36) 


(37) 


In  order  to  solve  (32)  and  (33)  with  the  appropriate  boundary  conditions, 
the  x-derivatives  are  approximated  according  to  equation  (5). 


Equations  (32)  and  (35)  are  linearized  with  respect  to  Cg  and  its  derivatives 
in  Newton-Raphson  fashion.  For  Tg  and  its  derivatives  the  last  iterates  are 
taken.  This  results  in  a  linear  two-point -boundary- value  problem  in  C^, 
which  is  solved  by  the  RH-method.  Subsequently,  the  same  quasi  Newton- 
Raphson  process  is  applied  to  equation  (33)  and  (36)  to  solve  Te  and  its 
derivatives.  To  deal  with  the  steep  gradients  near  the  wall,  a  suitable 
transformation  was  introduced  to  stretch  the  wall  region: 

n  =  (a+y/ymajc)1/n,  a  =  •001»  n  =  40  (37) 

If  only  rough  guesses  for  and  Tg  are  available  at  the  start  of  the  compu¬ 
tation,  the  iteration  process  of  both  coupled  boundary-value  problems  suffers 

from  an  instability  which  could  be  cured  by  relaxing  the  newest  value  of  T  : 

e 

T  =  <x  T  +  (1-a)  T  ,  (38) 

e-.ee, 
s+1  s  s+1 

in  which  s  is  the  number  of  iterations  and  Tg  is  the  result  of  the  quasi- 
Newt  on- Raphson  process. 

For  the  other  quantities,  the  results  could  be  accepted  for  each  iteration 
step  without  relaxation.  It  appeared  that ,  at  the  starting  station  only  the 
relaxation  parameter  a  had  to  satisfy:  .95  <  ol<  1. 

The  set  of  equations  (32)  through  (36)  is  solved  with  the  tuned  program  for 
the  cathode  wall  under  the  following  conditions:  current  crossing  the  chan¬ 
nel  Jy  =  -JO^A/m'  ,  the  emission  current  =  -1.03  10^A/m",  magnetic  field 
strength  B=1  Tesla.  It  is  assumed  that  the  gasdynamic  quantities  satisfy  the 
equation  of  state  and  at  the  centerline  satisfy  also  the  simplified  continuity 
and  momentum  equation:  p  =  yflRT, 

df  l/’cV  "  °’  (»> 

du 

P  u  T~~  +  -  J  .B. 

a  c  c  dx  dx  y 

T,  /O/yO, ,  u/u.,  are  given  functions  of  y  only  (fig.  4),  while  no  pressure 

variation  in  y-direction  is  present.  At  the  duct  entrance,  x=0,  the  following 

values  are  given:  velocity  at  the  centerline  UCQ=  400  m/s,  temperature  at  the 

centerline  T  =  1 890  K,  temperature  at  the  wall  T  =1600  K  and  pressure 
co  wo 

p  =1.5  Bar . 
o 

Fig.  5  shows  the  chosen  variation  of  u  and  P  along  the  duct.  At  the  initial 
position  x»0.1  the  starting  profiles  for  Tg  and  Cg  were  found  from  a  simpli¬ 
fied  form  of  the  continuity  and  energy  equation  (32)  and  (33),  in  which  con¬ 
vective  terms  arc  neglected. 
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In  figures  6  and  7,  results  are  depicted  of  C  —  and  T^-profiles.  In  figure  8 
a  detailed  picture  is  shown  of  the  effective  electrical  conductivity 
The  electron  concentration  shows  sharp  peaks  close  to  the  wall.  This  effect 
is  caused  by  the  emission  of  electrons  by  the  cathode  (Jem)*  If  should  be 
emphasised  that  the  model  does  not  permit  the  existence  of  space  charge, 
i.e.  the  number  density  of  ions  equals  the  number  density  of  electrons.  The 
electron  temperature  Te  increases  with  x,  due  to  the  decrease  of  (fig-  b) 
that  causes  less  energy  transfer  from  the  electrons  to  the  heavy  particles 
by  means  of  collisions.  This  effect  also  gives  rise  to  an  increase  in 
electron  density.  The  behaviour  of  T^  near  the  wall  is  less  obvious.  The 
model  of  the  boundary  conditions  at  the  wall  yields  for  the  energy  of 


electrons  passing  the  boundary  between  sheath  and  continuum  plasma: 

(Ee)  =  2kT  +eAjP  ,  whereas  the  mean  thermal  energy  of  electrons  in  the 
continuum  is  given  by  (Ee)  =  /2  kT  .  The  influx  of  electrons  from  the 
cathode  into  the  plasma  is  coupled  with  a  net  influx  of  energy  to  the  elec¬ 


tron  gas,  which  is  thermalized  by  electron-heavy  particles  collisions,  thus 
increasing  the  local  electron  temperature.  More  refined  models  for  the  plas¬ 
ma-electrode  interface  are  developed  (ref.  16,  17)  but  should  be  accompanied 


by  the  introduction  of  as  refined  and  more  complicated  models  for  the  sepe- 


rate  components  of  the  plasma.  The  computation  from  x=0.1  until  x=2.0, required 


190  sec.  computer  time  at  a  CDC  CYBER  72. 


5.  CONCLUSIONS 

Two  simple  methods  have  been  presented  for  the  solution  of  parabolic  diffe¬ 
rential  equations.  Both  methods  are  fourth-order  accurate  in  the  space-like 
dimension  and  second-order  accurate  in  the  time-like  dimension.  In  what  is 
called  the  L-method,  we  increased  the  accuracy  of  a  cubic-spline  method  by 
estimating  the  second-order  truncation  error  by  diff erentiating  the  basic 
equation  twice.  In  the  RH-method,  a  local  Richardson-extrapolation  is  em¬ 
ployed,  before  proceeding  to  the  next  x-station,  to  obtain  fourth-order 
accurate  values  at  every  two  points;  subsequent  quintic  Hennite  interpola¬ 
tion  is  then  used  to  increase  the  accuracy  of  the  remaining  points.  The 
boundary-layer  flow  near  a  stagnation  point  has  been  computed  with  both 
methods,  and,  for  comparison,  also  with  a  cubic-spline  method  and  with  a 
pure  difference  method.  Prom  a  von  Neumann  stability  analysis  applied  to 
the  internal  points,  it  is  inferred  that  the  RH-method  is  unconditionally 
stable.  The  D-method  has  the  same  stability  characteristics  as  the  second- 
order  cubic-3pline  method,  which  is  also  unconditionally  stable. 
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The  electron  gas  equations  in  an  MHD— channel  are  considered  to  be  of  suf¬ 
ficient  complexity  to  show  the  usefullness  of  the  RH— method  for  practical 
parabolic  problems.  The  equations  are  highly  non-linear,  and  the  boundary 
conditions  at  the  cathode  wall  are  functions  of  the  electron  gas  concen¬ 
tration  and  temperature  as  well  as  of  their  derivatives.  The  corresponding 
computer  program  required  some  tuning,  which  included  a  transformation  to 
stretch  the  y-coordinate  near  the  wail,  and  the  introduction  of  a  relaxa¬ 
tion  process  to  cope  with  the  high  non-linearity  of  the  enex-gy  equation.  The 
results  suggest  that  possibly  a  more  refined  model  for  the  plasma— electrode 
interface  is  to  be  considered  in  the  future.  The  RH-method  can  be  regarded 
as  a  suitable  tool  to  solve  boundary— layer  duct  flows  as  well  as  other 
complicated  parabolic  problems. 
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fig.  1  Schematic  view  of  an  MHD  generator  channel 


fig.  2  Accuracy  for  boundary  layer  flow; A  ,  difference  method;^,  cubic- 
splines;  0,  reduction  of  truncation  error  by  locally  applied 
Richardson-extrapolation ;  +, reduction  of  truncation  error  by  dif¬ 
ferentiation;  N*number  of  points  along  normal 


fig.  3  Eigenvalues  of  amplification  matrix  for  (k/h^) 6  =.l 


fig.  5  Variation  of  velocity  and  density  along  the  centerline  of  the 
duct  under  influence  of  the  Lorentz-force 


fig.  6  Electron  concentration,  (^-profiles  along  duct 
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ae  wave  number.  Substitution  of 
1  S'!-  and  5li_1  which 

l  by  the  introduction  of  5  2: 
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and  §  =(  §1,  £2)  , 


